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I. INTRODUCTION 



(N 

o 
o 

' We apply the synergetic elimination procedure for the stable modes in nonlinear delay systems 

}_( , close to a dynamical instability and derive the normal form for the delay-induced Hopf bifurcation 

■ in the Wright equation. The resulting periodic orbit is confirmed by numerical simulations. 

PACS numbers: 05.45.-a, 02.30.Ks 

in 
Q 

U, Within the last decades synergetics has provided powerful concepts and methods to describe selforganization 
P5 processes in various branches of science ^, ^, ^ |^. The spontaneous formation of spatial, temporal or functional 
' patterns in complex systems has been successfully investigated by working out general principles and by mapping 
them onto universal mathematical structures. The important result is due to the fact that in the vicinity of a 
dynamical instability the high-dimensional set of nonlinear evolution equations modelling a complex system on 
a microscopic or a mesoscopic scale can be approximately reduced to a low-dimensional set of order parameter 
^ ' equations describing the evolving pattern formation on a macroscopic scale. To obtain such a simplified, reduced 
10 • description of selforganization processes, the synergetic system analysis proceeds as follows. A linearization of the 
' evolution equations around a stationary solution shows that a dynamical instability is always accompanied by a 
time-scale hierarchy between numerous fast modes s and few slow modes u. A rigorous treatment of the full nonlinear 

■ evolution equations in the vicinity of the dynamical instability leads to a characteristic interdependence between 
[ both hierarchy levels which may be illustrated by a circular causality chain. On the one hand the slaving principle of 
' synergetics states that the numerous fast modes s quasi-instantaneously take values which are prescribed by the few 

slow modes u according to s{t) = h(u(t)) with the center manifold h(u). On the other hand an adiabatic elimination 
^ . of the fast enslaved modes s yields equations for the slow order parameters u which depend, in general, on the center 
^ manifold h(u) due to the nonlinear feedback. 

^ ' In its original formulation, the synergetic system analysis was developed for complex systems which can be modeled 
by ordinary and partial differential equations as well as their stochastic generalizations. Some time ago, the general 
concepts and methods of synergetics have been extended to delay differential equations to deal with dynamical 
I instabilities which are induced by the finite propagation time of signals in feedback loops |^ . Taking into account 
the infinite-dimensional character of a delay system |7| |^, the adiabatic elimination of the stable modes leads to a 
low-dimensional set of order parameter equations which turn out to be of the form of ordinary differential equations, 
i.e. they no longer contain memory effects. The predictions of the synergetic system analysis have been quantitatively 
tested by investigating the delay-induced Hopf bifurcation of the electronic system of a first-order phase-locked loop 
(PLL) [|6|. The periodic orbit which results from the corresponding order parameter equation near the bifurcation 
point has been confirmed by both a multiple scaling procedure and numerical simulations |l^ . Although this appli- 
cation exemplarily proves the order parameter concept for delay systems, it does not allow to draw conclusions about 
the slaving principle. As the lowest nonlinear term in the scalar delay differential equation of the PLL is a cubic one, 
the center manifold does not influence the order parameter equation of the Hopf bifurcation in the lowest order. In 
order to check for delay systems both ingredients of the circular causality chain, i.e. the order parameter concept and 
the slaving principle, it is thus indispensible to study a scalar delay differential equation with a quadratic nonlinearity. 
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A candidate for such a study is provided by 

d 



dt 



z{t)^R [z{t)-z{t-rf] . (1) 



With vanishing time delay r it represents a system which is named after the Belgian mathematician P.F. Verhulst 
from the 19th century It is used as a simplified model for the population dynamics of a species in an environment 
with limited food supply [|l2j. The synergetic system analysis for the Verhulst system with time delay (|^) has been 
already performed in Ref . There also the well-known equation of Wright 

^z(t) = -Rz{t-T)[l + z{t)] (2) 

has been treated, where R denotes a system parameter and r a delay time. This delay differential equation is mentioned 
by Wright jTslas arising in the application of probability methods to the theory of asymptotic prime number density. 
Cunningham ||14| depicts it as a "growth equation" representing a mathematical description of a fluctuating population 
of organisms under certain environmental conditions. In addition, it may describe the operation of a control system 
working with potentially explosive chemical reactions, and quite similar equations arise in economic studies of business 
cycles. Performing an appropriate scaling of time 

t = Tt', z'{t') = z{Tt') (3) 

converts the Wright equation (j^) to its standard form with the control parameter 

R' = tR. (4) 

Thus varying the delay time r corresponds to changing the control parameter R' . By omitting the prime ' for the 
respective quantities, the standard form of the Wright equation reads 

^z(i) = -Rz{t-l)[l + z{t)] . (5) 
In this paper we restrict ourselves to analyse this standard form of the Wright equation. 

The Wright equation (|^) shows a delay-induced instability, namely a Poincare-Andronov-Hopf bifurcation at the 
critical value 

Rc = ^ (6) 

of the control parameter R. In Ref. it is shown that the oscillatory solution in the vicinity of this instability, i.e. 
the emerging limit cycle, can be calculated approximately using the method of averaging. This approximation reads 
in the lowest order 



z{t) = A ^i? - I cos + O (i? - I) , (7) 
where the amplitude A has the value 



In Section ^ we start with a linear stability analysis of the Wright equation (j^) which confirms, of course, the delay- 
induced Poincare-Andronov-Hopf bifurcation when the control parameter R approaches the critical value (|^) . Near 
this instability we perform a nonlinear synergetic treatment in Section III and study in detail how the center manifold 
influences the order parameter equation. In Section IV the resulting order parameter equation is transformed to the 
normal form of a Hopf bifurcation, where the emerging periodic orbit is determined one order higher than the lowest- 
order result and (^). The numerical investigations of Section V confirm the emerging periodic orbit, furthermore 
we discuss the global bifurcation scenario of the Wright equation (j^). 
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II. LINEAR STABILITY ANALYSIS 



The solution of the delay differential equation (^) for times t >0 depends on the initial values of the function z(t) 
in the entire interval [—1,0]. In order to properly define such an initial value problem, Hale and Krasovskii 
proposed to transform the equation of motion for a function z{t) in the usual state space T to the extended state 
space C of continuous complex valued functions Zt(0) which are defined on the interval [—1,0]: 



(9) 



Following the notation of Ref. 
function z(t) € T through 



but also the nonlinear functional 



-1 < e < 
e = 

we introduced not only the new function zt (z C, which is connected to the original 



zt{e) ^ zit + e) , 



-i< e < 0. 



2 „o /.O k 

fc=i-^"i -^^1 1=1 



with the two scalar densities 



The stationary states of this system 



cj(2)(e,,e2) = -i?^(ei + 1)^(62). 



(10) 

(11) 

(12) 
(13) 

(14) 



are candidates for the reference state from which we start our further investigations. For the main body of the article 
we focus our attention on the stationary state z^^g^j and choose it as the reference state. The other stationary state 
zll^^ will be discussed in Section V together with the global bifurcation scenario of the Wright equation (||). 

Then we linearize the system (^ with respect to the stationary state zl^.^^. = by using the decomposition 

zt(e) = zLt + Ct(e), -i<e<o. (15) 



This leads to the following linearized equation of motion for the deviation Ct(Q) from the stationary state zl^. 



dt 



^C.(e) , -i<e<o 



where the linear functional is given by 



dew(e)Ct(e) 



e = 



with the scalar density 



Inserting the solution ansatz 



ST[zt] 



Sztie) 



-Rd{e + i) 



ztiey- 



Ct (6) =.0^(9)6^', 



4 < e < 



into (pj| ) leads to the eigenvalue problem of the infinitesimal generator Q^: 



= 0: 

(16) 

(17) 

(18) 

(19) 
(20) 
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Taking into account the definition of Ql in (|T6|), the eigenfunction (l)^{Q) is determined to be 

cj)^{Q)^ Nxe^'^ , -l<e<0, (21) 

and the eigenvalue A follows from 

A = L(A), (22) 

where L(A) is defined by 



L(A) / deuj{e)e^^ . (23) 



-1 



Using the scalar density (^^, we obtain the following transcendental characteristic equation: 

-Re-^-\ = 0. (24) 
Thus the spectrum of the linear operator Q l has the following properties Q : 

It consists of a countable infinite number of eigenvalues which cumulate for 5R(A) — oo. 
It is confined by an upper threshold for the real parts of the eigenvalues. 



At the bifurcation point, i.e. the instability, some of the eigenvalues reach the imaginary axes and thus become 
unstable. 

Further properties of the eigenvalues of the characteristic equation (^) follow from the Hayes theorem which can be 
found in Ref. Q . It states that all solutions of the transcendental equation 

p + qe~^ - X = (25) 



possess a negative real part if and only if (a) p < 1 and (b) p < —q < \/ a\ + p^. Here ai represents the solution 
of the transcendental equation oi = ptan(ai) which lies in the interval [0, tt). For the special case p = one can 
show that oi is equal to 7r/2. The shaded region in Fig. |l| represents that region of the parameter space p, q where 
both conditions of the Hayes theorem are fulfilled. The upper boundary line stems from (a) p < 1 and (bl) p < —q, 
whereas the lower boundary line follows from (a) p < 1 and (b2) —q < \/ a\ + p"^. 



Comparing (|24| ) with (g5|) we obtain the identification p = and q = —R. Changing the control parameter R 
from to 7r/2, the corresponding point in the parameter space p,q moves along the g-axis from the point q = to 
q = — Tr/2. At this critical value it reaches the boundary of the shaded stability region, i.e. no longer all solutions of 
the characteristic equation ( p^ ) have a negative real part. Therefore an instability occurs at Rc = i: 12. 

Fig. ^ confirms this result by illustrating the movement of the ten solutions of the characteristic equation ( p^ with 
the largest real part when the control parameter is increased from to 7r/2. The eigenvalues were obtained with a 
Newton algorithm and the control parameter R was increased in equidistant steps. For i? = there exists only one 
real eigenvalue as the linearized delay differential equation ( p^ degenerates to an ordinary differential equation. 
For i? > this eigenvalue remains real and becomes negative. Furthermore, a countable infinite number of conjugate 
complex eigenvalues and another real eigenvalue emerge from an infinite negative real part. At the value i? = 1/e 
both real eigenvalues meet at the point (—1/0). They are converted to a pair of conjugated complex eigenvalues for 
1/e < R < 7r/2. These two complex conjugated eigenvalues have zero real part at the instability Rc = 7r/2, which 
thus represents a Hopf bifurcation. We can further analyze this instability by introducing the smallness parameter 

^^R_Rc ^ R = Rc{l + £) (26) 

for the deviation from the critical control parameter Rc = Tr/2. In particular, we can determine both eigenvalues 
Aj(e) with nearly vanishing real part at e ~ from the characteristic equation (p^): 

AiW^^.±^iJ.(l + j^.)+0(.')^ (27) 
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FIG. 1: Within the shaded region both conditions of the Hayes theorem (a) p < 1 and (b) p < — g < ^/of^Tp^ with 
ai = ptan(ai) lying in the interval [0, tt) are fulfilled. Linearizing the Wright equation (^) around the stationary state z't^t — 
we have p = and q — —R. When the control parameter R is increased from to n/2, we move in the g,p-plane along the 
arrow within the shaded stability region. Thus all solutions of the characteristic equation have negative real part. The 
boundary of the shaded stability region is reached at the instability Rc = n/2. 




FIG. 2: The movement of the ten solutions of the characteristic equation ( p^ with the largest real part when the control 
parameter R is increased from to n/2. At the instability Rc = n/2 two complex conjugated eigenvalues have zero real part. 
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In the vicinity of the instabihty e « we read off from Fig. g that only the two eigenvalues ( |27| ) have nearly vanishing 
real part, all other eigenvalues have a large negative real part: 



3fi[A^(£wO)] wO; 5R[Ai(£wO)] <0, j = l,...,oo. 
This characteristic property of the linearized system ( p^ ) leads to a the time-scale hierarchy 

1 1 



5R[A±(e«0)] 



:$>n = 



5R[Ai(e«0)] 



J = 1, . 



(28) 



(29) 



Thus the infinite-dimensional extended state space C decomposes in an 2-dimensional subspace U of the linear unstable 
modes and a remaining infinite-dimensional subspace S of the linear stable modes As a consequence, the extended 
state function zt{Q) can be decomposed near the instability according to 



ztio) = zLt + + St (6) = u^io) + st{e) 



i< e < 0, 



(30) 



as we have zl^.^^ ~ 0. Here ^((9) and St(6) denote the respective contributions of Zt{Q) in the subspaces U and S. 
In order to project into these subspaces we need the linear unstable modes 



(6) = N.±e 



i< e < 



(31) 



of the system ( |16D which have been already determined in (|21|). However, their knowledge is not sufficient as the 
infinitesimal generator Ql is not selfadjoint. Therefore we need also the linear unstable modes 



-AS 



< s < 1 



of the adjoint system 



dt 



where the linear functional 



contains also the scalar density (p_8D. Indeed, the solution ansatz 




< s < 1 



converts ( |33|) to the eigenvalue problem 



AV't^(s) = (cjt ^tA) (s). 



< s < 1 , 



(32) 



(33) 



(34) 



(35) 



(36) 



which is adjoint to (20). Note that ip'''^{s) and Cj(s) are elements of the dual extended state space which consists 
of continuous complex valued functions on the interval [0, 1]. The relation between both extended state spaces C and 
C^^ is defined by the bilinear form 



=v;t (0)0(0) -y de j ds^Hs-e)uj{e)(t){s). 

-1 

Using this bilinear form one can show that the eigenfunctions (pl|) and (p2) are biorthonormal: 



(V^tA^I^A,, 



This determines the yet unknown normalization constants to be 

1 



^Aj = 



(37) 



(38) 



(39) 
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so they reduce near the instabihty because of (p7|) to 



N.± = , +0{e) . 



(40) 



Furthermore, the bihnear form (37) allows to define the projector into the 2-dimensional subspace U of the unstable 
modes: 



Zt 



(41) 



Correspondingly, the projector into the remaining infinite-dimensional subspace S of the stable modes reads 

Vs' = [J-Vu)' ■ (42) 
Applying the projector Vu to zt £ C leads to ut £U according to 

ut{e) = {Vuzt) (6) - J2 W-^"" ' (43) 

where the amplitudes of the linear unstable modes 0'^" (Q) are defined by 

^±(t)=(^tAj 



(44) 



Later on, these amplitudes represent the order parameters which indicate the emergence of an instability. Analogously, 
the projector (O) leads to the stable modes 



stie) ^ (Vszt) (e) , -i<e<o. 



(45) 



III. NONLINEAR SYNERGETIC ANALYSIS 



After having performed a linear stability analysis around the reference state zl^.^^. = in the vicinity of the instability 
Rc = 7r/2, we now return to our original nonlinear evolution equation (^) in the extended state space C. We proceed 
by decomposing the generator Q into its linear part Ql and a remaining effective nonlinear part: 



■yM^) = {gLZt){e) + Xoie)T^"[zt], -i<e<o. 

at 



Here we introduced the scalar function 



and the effective nonlinear functional 



Ao(e) 



, -1 < e < 

1 , e = 



^0 

dOi / de2u;^^\ei,e2) zt{Qi)zt{e2) 
-1 J-i 



(46) 



(47) 



(48) 



with the scalar density (|13|). Using the projectors ( [4l[ ) and ( [f2[ ) and their properties (|43|)~(}45|), we can investigate the 
respective contributions of the order parameters u^t) and the linear stable modes St G S to the nonlinear dynamics 
(Eq). Thus we obtain the following system of coupled nonlinear mode equations: 



3 = ± 



at 



J2^Kui{t) + st 
j=± 



(49) 
(50) 
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It is still exact and describes completely the nonlinear dynamics. However, a solution to these equations can be only 
found by means of an approximation method. Such a well-established approximative solution is provided by the 
slaving principle of synergetics ||, ||, ^, |5) . 



To this end we start with the time-scale hierarchy (29) near the instability which leads to the fact that the dynamics 
of the stable modes st & S evolves much faster than the order parameters (i) . In Rcf. it has been shown for 
a quite general class of delay differential equations that such a time-scale hierarchy leads to a slaving of the stable 
modes, i.e. the numerous fast modes St d S quasi-instantaneously take values which are prescribed by the few slow 
order parameters u^{t). In our context, the slaving principle states mathematically that the dynamics of the stables 
modes St € S is determined by the center manifold /i(0,u+,m~) according to 



stie) = h{e,u+{t),u-it)) . 

Inserting this ansatz in ( ^0|) leads to an implicit equation for the center manifold ft,(8,u+,' 



(51) 



dh{e,u+{t),u-{t)) 



J2 (i>^'-u\t) + h 

3 = ± 



■cff 



(52) 



It can be approximately solved in the vicinity of the instability as follows. We assume that the order parameters 



(53) 



' {t) possess a certain dependence on the smallness parameter (^ which is typical for a Hopf bifurctation 

u±(t) =0(ei/2). 

Furthermore, we perform for the center manifold h{Q,u^ the lowest-order ansatz 

jl=±j2 = ± 



(54) 



as r = 2 is the order of the effective nonlinear functional ( psf ) . From ( p3| ) and ( |5^ ) follows then in lowest order of e 
that the effective nonlinear functional JF*^^ in ( p2|) can be approximated by 



where the coefficients -Fj/J^ read 



■cff 



J2 (t)^'-u^{t) + h 
j=± 



ji=±i2=± 



(55) 



= j\ dQi y° dQ2 c^^'^ (Oi, 62) (ei)(/.^i' (62) . 
Taking into account (|l^) and (^l]), these coefficients turn out to be 



Ff, =F^1_* = -RNl+e^^i , Ff_ = F^% = -RN^ + N.-e'^i 



(56) 



(57) 



As a consequence, we conclude from ( p2| ) in lowest order of s that the coefficients Hj-^j^ (8) of the center manifold 
(U) are given by 



where the coefficients Kj-^j^ (Q) follow from 



Knj. (0) = {{Gl - I] ' {VuXo - Xo)) (6) 



(58) 



(59) 



with the abbreviation 



fe=i 



(60) 
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In Ref. 1^ it is shown that the operator [^^l — ^] ^ has the exphcit representation 

e / e 

{[Gl - 1]-' x) (e) = / ds p>(«-^)x(s) + [L{i) x(o) - fde f ds e^(e-^)c.(e)x(.s) 



(61) 



where L{X) is aheady defined in (^3|) . After some calculation, which also involves (^8|), (|T]), (^2|), (^ij), and (|47|), it 
thus follows that the coefficients (|9p are given by 



Ar2. 



Thus we obtain together with (|18D, (g3|), and (|60| 

K+-(e) = i^-+(e) = - 





iV2 e^"Q 


f — 


e2Aie 




A+ 


Xu 2Au 


Re- 


-2AJ + 2A+ 


)e 


Au 


A+ ' 


Re- 


(A++A-) ^ 


A+ + 



This completes the lowest-order result for the center manifold h{Q,u^,u ) which is given by (p^, (| 
and del). 



(62) 



(63) 
(64) 

del), 



Thus we can now consider the order parameter equation (^). In lowest order in s we take into account (|13[), (p3^), 
(^), (^, and ( ^ so that it reduces to 



dt 



=(t)=Aju±(i)i?iv,±n 



1=1 



where we set 



i=± ji=±i2=± 



-1 , ; = i 

,1 = 2 



(65) 



(66) 



Note that the order parameter equation ( pSj ) turns out to be an ordinary differential equation, i.e. it no longer 
contains memory effects. Furthermore, we observe that the center manifold explicitly enters the order parameter 
equation ( |65|) as a direct consequence of the quadratic nonlinearity of the Wright equation (H) . We remark that this 
effect, which is essential for the present synergetic analysis, was neglected in the neurophysiological study in Ref. p9[ . 
In the subsequent section we show how the order parameter equation ( |65| ) is converted to the normal form of a Hopf 
bifurcation. 



IV. NORMAL FORM 



Now we perform a nonlinear transformation of the order parameters which eliminates those terms which are irrel- 
evant for the normal form of a Hopf bifurcation. As far as the so-called near identity transformation and the theory 
of normal forms in general is concerned, we refer to the Refs. The terms in (|65| ) which are relevant for 
the normal form of a Hopf bifurcation read 



j^u^it) = Xtu^{t) + qtu^{tf + qtu^{t)u^{t) + qtu^itf + k^u^{tfu^{t) 



(67) 
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as we can neglect quartic terms and nonresonant cubic terms due to the rotating wave approximation. The respective 
coefficients in Eq. (|67|) are given by 



1^ 



1i 

it 
kt 



-i?7V^±0^J (-1)0^^(0) 



-Aj 



-i?A^^±(/)^J (-1)0^J(O) = -RN^±N^T^e-^i 



{(/.^J (-1) [H+-{0) + i/-+(0)] + {-1)H++{Q) + 
cp'iiO) [i/+_(-l) + i/_+(-l)]+0^J(O)i/++(-l)} , 



(68) 
(69) 
(70) 

(71) 



where we did not write down the exphcit form of kf^ for simphcity. Then the previous order parameters u^{t) are 
transformed to new order parameters v^{t) by the near identity transformation 



±„.±i 



(72) 



with the yet-unknown coefficients aj, a^, and at- As the u^(t) are small quantities in the vicinity of the instability, 
the same holds for the v^{t). Inserting ( |72|) in (|67|), we obtain a system of ordinary differential equations of the form 



M{t) ■ 



d fv+it) 



w+it) 
w-{t) 



dt yv'it) 

where the matrix M{t) is defined by 

_ f 1 + 2a^v+ {t) + a+ir (t) a+v+ [t) + 2a+u" (<) 



M(t) 



^1 ir (<) + 2a2 v+ {t) 1 + 2^0 {t) + v+ {t) 



(73) 



(74) 



For simplicity, we don't write down the explicit form of w^{t) and w (t), but we note that they contain v^{t) and 
v~'{t) at least in first order. Thus we obtain from (73) 



d_ /v+{t) 
dt [v-it) 



M-\t) 



w+it) 
w-{t) 



with the inverse matrix 



M{ty 



1 



M22{t) -Mi2{t) 



DetM(t) V ~M2i{t) Mn{t) 
where the determinant has the form 

Det M{t) = 1 + v+{t) (2(1+ + a~) + 2v+{t)v~{t) {a+a~ - a^a~) + 2v+[tf {a+a~ - a+a") + c.c. 
Expanding the right-hand side of ( ^5| ) in powers of v+{t) and v~ (t) up to the third order, we yield 
d 



(75) 



(76) 



(77) 



dt 



v^{t) = A±w±(t) + {qf - a±A±) v^ity + {qf - afX+) v+{t)v-{t) + qf + a±(A± - 2A+) 



v+ity 



kf + qoiaf - aj) - q^af + qfaf + 2qfa+ - 2q+af + afafX^ + 2af a^(2A± - A+) 
+a^at{X^ +2X+)\v^{tfv+{t). 



(78) 



Now we can fix the yet- unknown coefficients , af^ , and of the near identity transformation ( fz^ ) in such a way 
that all quadratic terms vanish. This leads to the conditions 



an 



4 



Xu Xu 2Xu — Xu 

Thus Eq. (|7^) reduces to the normal form of a Hopf bifurcation 

|«±(^) = A^^;±(^) + 6±«±(^)V(^), 



12 



(79) 



(80) 



11 



where the Hopf parameter is given by 

kt + 



An Atf (2Au — Am ) 



Taking into account (|2|), @, (|T|), and @ as well as @, (|8|), (|3|), and (||) together with 
parameter reads in the vicinity of the instability: 



Rc 



[{SRc - 1) ± liRc + 'S)]+0 (e) 



5(l + i?2)f 
Performing the ansatz 

the normal form (|8^) is transformed to polar coordinates 

d 



dt 
d 

di 



r{t) = r{t)[^{\^)+^{b^)r{tf 



Thus taking into account (E7f) and (p3) near the instability, we result in the oscillatory solution 



^stat 



SR(A^) 
" 5R(5±) 



5i?c 



= ±[Q(A^)+3(6±)rf 



3i?^ - 1 



(81) 

ID, this Hopf 
(82) 

(83) 

(84) 
(85) 

(86) 
(87) 



In order to compare this result with numerical simulations, we have to convert this oscillatory solution back to the 
original state space F. At first we observe that we obtain for e C from (^o|), (^), (51), and (|j) near the instability: 

i=± ji=±j2=± 

Taking into account the near identity transformation (^) together with (^), this yields up to the first order in e: 



zt{e) = N^+e^i^v+{t) + iV^-e^" ^^-(t) + ao(e)?;+(0 + ai{e)v+ it)v- (t) + a2ie)v-{t) 
where the coefficients ao(6), ai(9), 02(6) read 

ao(e) = \N^+e^^^a+ + N^-e^^^a^ + H++{e) 

02(6) = [iV;,+ e^^®a+ +7V^-e-^"®ao +i/__(e) 

Due to the relation (^0|) between z{t) gT and zt € C, we conclude from (|8^): 

z{t) = A^;,+ w+(t) + ^A-«" W + ao(0)w+(0^ + ai(0)v+(t)w-(t) + a2{0)v-{tf 
Near the instability we obtain from (Bc 



1 



with some phase i/'i; whereas (B3)-(p5h leads to 



(89) 

(90) 
(91) 
(92) 

(93) 

(94) 

(95) 
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with the radius (^6|) and the phase 
Here the frequency turns out to be 



(p{t) = n{e)t + ipo ■ 



Rc 



3i?c- 1 

Furthermore, we yield from (p0[)-(p2[) by taking into account (|57|), (| 

1 



ao(0) 



ai(0) = 0, a2(0) 



(96) 
(97) 

33|), (p^), and ( p4[ ) in the lowest order of e 



V5(i + Rl) 



V5(i + Rl) 

where '02 denotes some phase. Thus we obtain the following result for z{t) G F near the instability 

z{t) = co(£) + ci(£) COS [ip{t) + V'l] + C2{e) COS [2ip{t) +^2]+0 (e^'j , 
where the respective coeffcients read 



Co(£)=0, Ci(£) = 2 



bRc r / ^ o ^R^ 



iRr 



3Rr 



1 



(99) 



(100) 



Now we compare the oscillatory solution (^ , (|^) , which was obtained from using the method of averaging, with ours 
(p6|), (p7|), (p9|), ( |100| ) by taking into account the critical value (^ of the control parameter. We conclude that both 
results coincide in the lowest-order £^/^, but our result is even correct up to the order e. 

From the near identity transformation (|7|) as well as from (|86| ) and (|9^) we conclude that the order parameters 
u^{t) turn out to be of the order £^/^. This result is consistent with our original assumption (|5^) which was the basis 



of our approximate solution of the implicit equation for the center manifold ( |52| ) in the vicinity of the instability. 
Thus our synergetic system analysis is justified a posteriori by self consistency. 

Note that the same perturbative result (p6|), (|97), (p9|), (IOC) for the oscillatory solution above the Hopf bifurcation 
can be derived with the multiple scaling method 1 9] . It represents a technical procedure to deduce the normal form, 
once the bifurcation type is known, by using the knowledge how the respective quantities depend on the smallness 
parameter e = {R — Rc)/Rc- Although the multiple scaling method has been originally developed for ordinary 
differential equations ||2^, |l, 
inRef. l2l). 



22 , it can be also applied to delay differential equations (see, for instance, the treatment 



V. NUMERICAL INVESTIGATION 



In order to numerically verify our analytical result, we integrated the underlying delay differential equation of 
Wright (|^). By doing so, we varied the control parameter R in the vicinity of the instability Rc = 7r/2 in such a 
way that the smallness parameter e — {R ~ Rc)/Rc took 200 equidistant values between 10~^ and 10~^. We used 
a Runge-Kutta-Verner method of the IMSL library as an integration routine and performed a linear interpolation 
between the respective values. In particular in the immediate vicinity of the instability the phenomenon of critical 
slowing down led to a transient behavior. To exclude this, we iterated the discretized delay differential equation for 
each value of the control parameter at least 10^ times. Afterwards we calculated the power spectrum with a complex 
FFT so that the basic frequency fl of the oscillatory solution could be determined with high resolution. Then we 
performed a real FFT with the period T = 27r/f2 of the simulated periodic signal z{t) — z{t + T): 



z{t) = ^+^[ak cos (kilt) + bk sin {kilt)] . (101) 



a-o 

k=l 

The Fourier coefficients follow from integrations with respect to one period T — 27r/51: 



oo oo 

ak = ^ J dtf{t) COS (knt) , fc = 0,l,... ,oo; bk = ^ J dt f (t) sin {knt) , fc = l,...,oo. (102) 
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FIG. 3: Frequency and Fourier coefficients cq, ci, C2 of the osciilatory soiution of the Wright equation after the Hopf 
bifurcation versus the smallness parameter e — {R~ Rc)/Rc- The interval of the smallness parameter [10~^, 10~^] was divided 
into 200 equidistant parts. 



From (101) follows then the spectral representation 

oo 

z{t) = Co + ^Ck COS {km + (l)k) (103) 



with the quantities 



k=l 



Co = ^ , Cfc = ^Jal + bl , 0fc = -artan — , k = 1, . . . ,00 . (104) 



Thus our analytical result (j^m, ( p9[ ) can be interpreted as the first terms within a spectral re pres entation (103), 
where the frequency = 2tF/T and the Fourier coefficients cq, ci, C2 are given by (p7| ) and ( 100|) . Numerically 
analysing the Hopf bifurcation with the FFT, the results for fl, cq, ci, C2 are plotted in Fig. ^ versus the smallness 
parameter e. Comparing the respective numerical and analytical results, we observe some deviations for small and 
for large values of the smallness parameter e. The former are due to the phenomenon of critical slowing down, 
i.e. the system stays longer in the transient state when the instability is approached, and the latter arise from the 
neglected higher order corrections in the analytical approach. Therefore we restricted our numerical analysis to the 
intermediate interval [10~^, 10~^] of the smallness parameter e. In Tab. | we see that the analytical and numerical 
determined quantities agree quantitatively very well. Thus our synergetic system analysis for the delay-induced Hopf 
bifurcation in the Wright equation is numerically verified. 

For the sake of completeness we have also investigated oscillatory solutions for values of the control parameter R 
which are larger than the critical one Rc = tt/2. Figure ^ shows that all these periodic solutions oscillate around 
the stationary state zl^^^^ — which becomes unstable at Rc — Tr/2. It turns out that a global bifurcation occurs for 
R^ = 3.247 as then the oscillatory solution comes close to the other stationary state zH^^^ = — 1 which turns out to be 
linear imstable for all values of the control parameter R > 0. Indeed, performing a linear stability analysis according 
to Section || around the stationary state zH^^ = — 1 leads to the characteristic equation 

R-X = 0, (105) 

so we have from Eq. (24) the identification p = R and q = (compare the shaded stability region in Fig. 111). 



VI. SUMMARY AND OUTLOOK 



In this article a linear stability analysis of the Wright equation (^) around the stationary state zl^^^ = showed 
that a delay- induced Hopf bifurcation occurs at the critical value Rc = tt/2 of the control parameter R. Within a 
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investigated 
quantity 



analytical 
expression 
axis intercept slope 



analytical 
value 

axis intercept slope 



numerical 
value 

axis intercept slope 



Rc 



Rc 



3Rc - 1 



1.5708 



^0.4231 



1.5707 



-0.4024 



co(e) 



0.0 



0.0 



-2 ■ 10" 



4- 10" 



ci(e) 



5Rc 



- 1 



2.9090 



0.5 



2.890 



0.4999 



C2(e) 



V5R, 
SRc-1 



0.0 



1.8923 



2 ■ 10" 



1.832 



TABLE L Comparison between the analytical and numerical values for the frequency Q and the Fourier coefficients co, ci, C2 
of the oscillatory solution of the Wright equation after the Ifopf bifurcation. 




FIG. 4: Oscillatory solutions of the Wright equation for three values of the control parameter R: 1.7 (dotted line), 2.0 
(dashed line), and 3.247 (solid line). 



subsequent nonlinear synergetic analysis we adiabatically eliminated the stable modes and derived the normal form 
of this Hopf bifurcation. It is explicitly influenced by the center manifold in the lowest order, as the Wright equation 
(^ has a quadratic nonlinearity. Solving the normal form we obtained a periodic solution above the Hopf bifurcation 
which was numerically verified. 

In contrast to the corresponding analysis of the electronic system of a first-order phase-locked loop with time delay 
this paper confirms not only the order parameter concept for delay systems. It also represents a successful test 
for the slaving principle of synergetics, i.e. for the influence of the center manifold on the order parameter equations. 
Thus the validity of the circular causality chain of synergetics (see Fig. ||) has been demonstrated for the Hopf 
bifurcation of a delay differential equation. 

It remains to investigate the circular causality chain also for other bifurcations. The Floquet theory for delay 
differential equations and thus the linear stability analysis for a periodic reference state was already established in the 
Rcfs. 24, 251]. However, a corresponding synergetic system analysis is still missing which derives the order parameter 
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few mode amplitudes 




i 


of the slow 






linear unstable modes 





order parameter u'^lt) 
center manifold 

5,(6) = h(e,u+(t),u-{t)) 



1 



order parameter equations 
Vjt^J (0) J^efF 0A,t„±(j) + J 



many fast 
linear stable modes s/(©) 



FIG. 5: Circular causality chain of synergetics for the Hopf bifurcation of a delay differential equation. On the one hand the 
center manifold of the slaving principle guarantees that many fast linear stable modes St (O) quasi-instantaneously take values 
which are prescribed by the few slow linear unstable modes M^(t). On the other hand the adiabatic elimination of the fast 
enslaved modes St(0) influences the resulting order parameter equation. 



equations and the normal forms for bifurcations of oscillatory solutions 
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